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ABSTRACT 

It appears that most stars are born in clusters, and that at birth most stars have 
circumstellar discs which are comparable in size to the separations between the stars. 
Interactions between neighbouring stars and discs are therefore likely to play a key 
role in determining disc lifetimes, stellar masses, and the separations and eccentricities 
of binary orbits. Such interactions may also cause fragmentation of the discs, thereby 
triggering the formation of additional stars. 

We have carried out a series of simulations of disc-star interactions using an SPH 
code which treats self-gravity, hydrodynamic and viscous forces. We find that interac- 
tions between discs and stars provide a mechanism for removing energy from, or adding 
energy to, the orbits of the stars, and for truncating the discs. However, capture during 
such encounters is unlikely to be an important binary formation mechanism. 

A more significant consequence of such encounters is that they can trigger frag- 
mentation of the disc, via tidally and compressionally induced gravitational instabil- 
ities, leading to the formation of additional stars. When the disc-spins and stellar 
orbits are randomly oriented, encounters lead to the formation of new companions to 
the original star in 20% of encounters. If most encounters are prograde and copla- 
nar, as suggested by simulations of dynamically-triggered star formation, then new 
companions are formed in approximately 50% of encounters. 

Key words: stars: formation - binaries: general - accretion, accretion discs - meth- 
ods: numerical - hydrodynamics - instabilities 



1 INTRODUCTION 

In relatively quiescent molecular clouds, the star formation 
process may be initiated quasistatically, by ambipolar diffu- 
sion; individual protostars then condense out indepedently 
of one another, via inside-out collapse. Star formation may 
also be initiated dynamically, and in this case multiple pro- 
tostars are formed; here we have in mind impulsive trig- 
gers such as c ollisions between clumps in molecular clouds 
( Pringle 1989 ) or the compression of material by expanding 
nebula, i.e. HII regions, stellar- win d bub bles, and supernova 
r emna nts (e.g. Elmegreen & Lada (L977)). Whitworth et al. 
( 1996 ) have argued that dynamically triggered star forma- 
tion is significantly more prolific than quiescent star forma- 
tion, and that it is much more likely to produce star clusters 
and binary systems. 

In quiescent star formation, a core within a molecular 
cloud collapses to a protostellar disc. Viscous forces within 
the disc then remove angular momentum from the disc mate- 
rial, allowing it to accrete onto a central protostar. This ac- 
cretion process continues until, by the time the protostar has 
reached the T-Tauri stage, most of the mass of the system is 



in the central star, with only a small fraction remaining in a 
vestigial disc. Discs are probably present around protostars 
for a large part of their pre-main-sequence lifetimes. 

In contrast, dynamically-triggered star formation is 
more chaotic. The simulations of dy nami cally triggered star 
for mation reported by Turner et al. ( 1995) and Whitworth et 
al. ( |199£| ) suggest that small- A'^ clusters of protostellar discs 
are formed. Impulsive interactions between these discs then 
occur on a timescale comparable with - or possibly even 
shorter than - the timescale for the disc to accrete onto the 
central star. The main purpose of this paper (Paper I) and 
its two companion papers is to explore the consequences of 
these impulsive interactions. 

Observations indicate that most stars are born in clus- 



ters (e.g. Lada et al. (1991)). Observations also suggest that 



most stars condense out of protostellar discs having radii in 
the range 100-lOOOAU, as compared to separations between 
stars in clusters of approximately 5000-lOOOOAU (Strom et 
al. 1989b, Lada et al. 1991). Therefore it seems unavoid- 
able that most protostellar discs will, during their lifetime, 
undergo at least one encounter with a neighbouring proto- 
stellar disc or star. 
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It has been proposed (Larson 199C) that in dense stellar 



clusters, capture due to gravitational drag caused by star- 
disc interactions could be an important binary formation 
mechanism. However, this mechanism cannot explain the 
formation of binaries wi th per iastra 1z 1 AU or ^ 1000 AU; 
and Clarke fc Pringle (1991a) have shown that, except in 



of capture as a binary formation mechanism. He finds that 
encounters within 2 disc radii can lead to the stripping of 
approximately half of the disc mass, but that capture at this 
epoch cannot be a dominant binary formation mechanism. 
Simulations of the effect of a perturber on a disc have 



also been carried out using finite difference methods (Ko- 



the dersest stellar environments, capture cannot lead to a rycansky & Papaloizou 1995). These simulations only treat 



significant binary fraction. 



Clarke & Pringle (1991b) and McDonald & Clarke 



(1995) have modelled the evolution of an initially-bound 
small-Ai' cluster of stars, treating the dissipative effects of 
circumstellar discs with a simple prescription for the energy 
loss during an encounter. They find that the discs have a 
significant influence on the evolution of the system, leading 
to a considerable number of binary and multiple systems. 
However, their results are sensitive to the prescription for 
interactions used. 

For weak encounters, in which the periastron distance is 
much greater than the disc radius, it is possible to use linear 
perturbation analysis to calculate the energy and angular 
momentum transfer between the disc and the orbit of the 



two stars (Ostriker 1994). At these large separations energy 



and angular momentum are in general transferred from the 
disc to the stars. Thus long-range encounters do not provide 
a mechanism for capture, but can give an enhanced accretion 
rate within the disc. 

Encounters in which the perturber passes close to, or 
penetrates, the disc are highly non-linear and cannot be fol- 
lowed analytically. It is therefore necessary to simulate such 
encounters numerically. Below we review the results of such 
simulations. 

Clarke & Pringle ( 1993) have used a reduced three-body 



method to investigate the response of an accretion disc to a 
perturbing star. Their discs were non-self-gravitating, pres- 
sure forces were ignored, and viscous forces were treated us- 
ing a pseudo- viscosity chosen to prevent fluid elements from 
passing through one another. They investigated three differ- 
ent cases. They found that for prograde, co-planar encoun- 
ters, the disc is tidally stripped to about half of the perias- 
tron distance, and the perturbing star captures a significant 
amount of the disrupted material. In a retrograde, co-planar 
encounter, the disc is disrupted to about the periastron ra- 
dius, and very little material is captured by the perturbing 
star. For orthogonal encounters, the disc is disrupted out- 
side of periastron, and the remaining disc is twisted. The 
perturbing star captures almost no material. 

A more detailed series of simulations has been carried 



out by Heller (1993), using smoothed particle hydrodynam- 



ics (SPH). The primary aim of this work was to calculate the 
tilting of the disc due to the passage of stars on non-coplanar 
orbits. He included pressure forces, using an adiabatic equa- 
tion of state for the disc. He investigated a larger parameter 
space than Clarke & Pringle, with periastron distances from 
0.5 to 2 disc radii, and a range of orbital inclinations, and 
allowed for the effect of the disc upon the orbit of the stars. 
He confined his simulations to a IMq star, surrounded by a 
O.IMq disc of radius 100 AU. This corresponds to an evolved 
system, rather th an the massive, extended discs suggested 
by Lin & P ringle ( |l99C| ). 

Heller (1995) has used the same method, this time with 
a locally isothermal equation of state, to investigate the ef- 
fect of encounters upon the disc material, and the possibility 



the case of long-range, co-planar prograde encounters. They 
find that the encounter removes angular momentum from 
the disc, leading to an enhanced accretion rate, and that 
energy is transferred from the dis c to t he stars - in agree- 
ment with the results of Ostriker (1994). 

In order to discover which parts of the disc are most re- 
sponsible for the transfer of en ergy and angular momentum. 
Hall, Clarke & Pringle ( 1995 ) have carried out a series of 
reduced 3-body simulations, in which hydrodynamic forces 
are ignored. They find that the transfers of both energy and 
angular momentum are usually dominated by the material 
that becomes unbound during the encounter. They also find 
that for prograde encounters, a corotation resonance causes 
energy and angular momentum to be transferred to the bi- 
nary orbit, so that for encounters with the periastron radius 
greater than about twice the disc radius, the binary may 
become less strongly bound, or even unbound. 

This paper is the first of three which describe simula- 
tions of encounters involving massive, extended protostellar 
discs. In this paper we describe our numerical method, dis- 
cuss its limitations, and present simulations of interactions 
between a discless star and a star possessing a massive, ex- 
tended disc. The second and third papers (Watkins et al. 
1997b,c) investigate interactions between two protostellar 
discs. 



2 NUMERICAL METHOD 

The method we use in our simulations is smoothed parti- 
cle hydrodynamics (SPH ) combined with treecode gravity 
(Hernquist & Katz 1989). The discs in the simulations are 



self-gravitating, and experience pressure and viscous forces, 
as well as the gravitational attraction of the stars. The stars 
are represented by particles that experience only gravita- 
tional forces, and which move in the combined potential of 
each other and the discs. Our code uses an adaptive smooth- 
ing length h, so that each particle has an approximately con- 
stant number (~ 50) of neighbours. For instance, the par- 
ticles in the initial disc have h ~ 80AU. Gravity is kernel- 
softened, so h is also the gravity-softening length for the 
particles. The stars have a constant gravity softening length 
of 50 AU. 



SPH is normally used with an artificial viscosity (Mon 



aghan 198!-), in which the bulk viscosity and shear viscosity 
are coupled together. A high bulk viscosity is necessary in 
order to prevent unphysical penetration of colliding particle 
streams, and hence to capture shocks, but this means that 
the resulting shear viscosity is too strong. 

Here we u se a prescription for the viscosity due to 
Watkins et al. (1996). This prescription allows the bulk and 



shear viscosities to be regulated independently, so that it 
is possible to have a bulk viscosity high enough to model 
the shocks in disc-star interactions, without having an un- 
physically high shear viscosity. We use a bulk viscosity 
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property 



value 



disc radius 
primary disc mass 
primary star mass 
secondary star mass 
disc density profile 
eccentricity of orbit 
initial separation 
shear viscosity 



lOOOAU 

O.5M0 

O.SMq 

1.OM0 

,.-3/2 

1.0 

5000AU 

Us = Csh/800 



Table 1. Physical parameters used in simulations 

Hb = Cshp/8, where Cs is the sound speed, h is the SPH 
smoothing length and p is the density. This viscosity is 
equivalent to the bulk compon ent of the g-t erm of the SPH 
artificial viscosity with a ~ 1 (Murray 1996). For the shear 
viscosity we use /is = Cshp/800, w hich i s equivalent to the 



Q- viscosity of Shakura & Sunyaev (197S) with ass ~ 10 



We also include the 0- term from the artificial viscosity with 
P ~ 1 (seeMonaghan (198£)). We have performed a variety 
of tests on this vis cosity prescription, and they are reported 
in Watkins et al. ( |l99e| ). 

Initially the numerical Reynolds number in the discs 
is ~ 3000(r/1000AU)"^/^. There is no indication of turbu- 
lence developing, but this may simply be a consequence of 
the low resolution. Isolated discs evolve secularly on a vis- 
cous times cale. and in the sense predicted by Lynden-Bell 
& Pringle (|l974|). 



2.1 Limitations of the method 

The disc is modelled using a single layer of 2000 particles, 
and so it is not resolved vertically. This means that for 
non-coplanar encounters, the propagation of sound waves 
through the disc may be inhibited. However, the interac- 
tions we simulate are dominated by gravitational torques 
and shock compression, so sound waves have insufficient 
time to have a significant effect on the outcome. We infer 
that the thinness of the discs does not corrupt the results 
obtained, but this can only be confirmed by carrying out 
simulations using a much larger number of particles, so as 
to resolve the discs vertically. 

We were obliged to use a relatively small number of par- 
ticles for three reasons. First, we needed to explore a rather 
large parameter space, and so the individual simulations had 
to be relatively short. Second, we needed to follow the con- 
densation of individual protostellar discs, and this incurred a 
rather short time-step particularly in the later stages of the 
simulations. Third, we have been unable to marry our new 
viscosity with our multiple time-step scheme, and so the sim- 
ulations become particularly compute-intensive as the new 
protostars condense out in the later stages. 

In the central regions of the disc, where the particles 
are orbiting very fast, particles may find neighbours with 
supersonic relative velocities. Indeed, they may even find 
neighbours on the opposite side of the disc. This leads to 
the particles in the central region experiencing an unphysical 
shear that causes them to decelerate and spiral in towards 
the central star much faster than they should do. This has 
the effect of altering the initial surface-density profile in the 
immediate vicinity of the star. In combination with grav- 
ity softening, the high shear also causes the central region 



to form a pseudo-lattice in solid-body rotation, as can be 
seen in Figure 3. This problem remains isolated to within 
~ 2h of the centre of the disc, and does not propagate out- 
wards. As a consequence, the inner 10% of the disc is not 
correctly modelled, and any fragmentation that might occur 
in this region is artificially suppressed. In order to alleviate 
this problem, it would be necessary to use a larger number 
of particles, so that h could be smaller and hence fewer par- 
ticles would have supersonic neighbours. The problem could 
also be overcome by the use of sink particles in the centre 
of the disc (Bate, Bonnell & Price 1995). 



3 PHYSICAL AND COMPUTATIONAL 
MODEL 

During the Class and Class I stages, young protostars 
are heavily obscured by their surrounding envelopes and are 
thus extremely difficult to detect. The vast majority of ob- 
servations of pre-main sequence stars are therefore of the 
later Class II and HI stages, i.e. T-Tauri stars. At this stage 
of the evolution, the central star has accreted most of its 
surrounding disc, so that the discs are typically of radius 
lOOAU and mass 0.001 to O.OIMq (e.g. Beckwith & Sar- 
gent 1993). However, observations indicate that by the time 
stars reach the T-Tauri stage, they already have a binarit y 
at least equal to that of main sequence stars (Mathieu 1994). 



Whatever determines the binary statistics of protostars must 
therefore occur at an earlier stage than the T-Tauri phase. 
For this reason, we have chosen to model the case in which 
the star and disc masses are equal, and the disc is extended 
(lOOOAU). This corresponds to a Class or young Class I ob- 
ject. The primary star and its disc are both of mass 0.5Mq, 
whilst for the star-disc encounters, the perturbing star has 
mass I.OMq, equal to the total mass of the system it is en- 
countering. The physical parameters used in this paper are 
summarised in Table [l|. 

Lin & Pringle (199C) found in their collapse calcula- 
tions that when a disc formed, its surface density profile 



began as E oc 



and evolved to E oc r 



-3/2 



This evolu- 



tion was driven by transport of angular momentum due to 
gravitational instabilities in the disc. The simulations pre- 
sented here therefore use an initial surface density profile 
E oc (r -I- ro)"^''^, where ro = 0.02rdisc is a small term to 
smooth out the central singularity. In the standard simu- 
lations, the disc is evolved in isolation before the start of 
the encounter, in order to remove any initial transients, and 
to allow the disc to reach a state in which it evolves on a 
secular, rather than dynamic, timescale. 

The discs are taken to have an isothermal equation of 
state at lOK. As the discs are optically thin to their mm 
cooling radiation, this should be a good approximation ex- 
cept in the very central regions, where radiation released 
from Kelvin-Helmholtz contraction of the central protostar 
will heat the disc. At densities of greater than lO"'^^ g cm"'', 
condensation is slowed down by making the equation of state 
polytropic with exponent 1.7. The results of the simulations 
can be rescaled to different physical units. For example, they 
are also applicable to 300AU radius discs which have a tem- 
perature of 33K. 

The orbit of the perturbing star is expected to be 
mildly hyperbolic (Larson 1990). There is very little dif- 
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Run 


"^'peri 1 disc 





dsOl-04 


0.5, 1.0, 1.5, 2.0 





ds05-08 


0.5, 1.0, 1.5, 2.0 


7r/4 


ds09-12 


0.5, 1.0, 1.5, 2.0 


7r/2 


dsl3-16 


0.5, 1.0, 1.5, 2.0 


37r/4 


dsl7-20 


0.5, 1.0, 1.5, 2.0 




Table 2. 


List of simulations 





ference energetically between mildly hyperbolic encounters 
and parabolic ones, so for simplicity all of the encounters 
modelled here are parabolic. We have carried out a suite of 
standard simulations varying the periastron distance, r^^ri^ 
and the orbital inclination, The simulations carried out 
are listed in Table ^. We begin the simulations with the 
perturbing star at five times the disc radius. 

We have also carried out a number of tests using higher 
resolution, starting the discs further apart, using different 
initial surface-density profiles, and omitting to settle the 
discs in isolation before the start of the encounter. No signif- 
icant differences in the results are seen. Appendix B presents 
results obtained when the disc is still modelled with a sin- 
gle layer of particles, but with more particles (11299 instead 
of 2000) and hence greater resolution. Appendix C presents 
results obtained when the disc is resolved in the third di- 
mension (parallel to the spin axis) using 13713 particles. In 
both cases the results are remarkably similar to those ob- 
tained with the standard resolution, and this gives us confi- 
dence that the basic physics is not being corrupted by lack 
of resolution. 

An additional test simulation has been carried out to 
investigate the secular stability of an isolated disc, and to 
determine whether the growth of intrinsic gravitational in- 
stabilities within the disc leads to its fragmentation. The 
evolution of the disc has been followed for 20 orbital periods 
at the outer edge. One- and two-armed spiral instabilities 
develop within the disc but act only to redistribute angu- 
lar momentum within the disc, leading to accretion onto 
the centre and the spreading of the outer regions of the disc. 
This is as expected, as the disc that we have chosen t o model 
has a m inimum value of the Toomre Q-parameter (Toomre 
1964 |) c| f 2.4, ai id instability to fragmentat ion is expected 
only for Q J: 1 (Laughlin & Rozyczka 1996). The relatively 



high Q of the disc means that we should not expect the disc 
to retain the r"^''^ surface density profile derived by Lin & 



Pringle (1990) 



4 SIMULATIONS 

The figures presented in this section show the results of two 
of the simulations conducted. The figures are grey-scale plots 
of column-density, in which the shading is logarithmically- 
scaled. These plots are overlaid with contours of constant 
surface density, that are equally separated in log-space. Oc- 
casionally an extra contour has been added in order to high- 
light an important feature of the figure, such as a shock or 
spiral arm. The contour levels used are given in the caption 



to the figure. The time at each figure, in units of lO** years, 
is also given, and is calculated from the start of the simula- 
tion, when the two stars are at a separation of 'hrdiac- The 
individual frames within a figure, when referred to within 
the text or a caption, are labelled from top left to bottom 
right in rows, so that frames 1 and 2 are on the top row, 
frames 3 and 4 on the next row and so on. 



4.1 run ds02 : prograde coplanar encounter 

Encounters with </> = are prograde, coplanar encounters. 
The results for run ds02 are shown in Fig. |l| The periastron 
distance is lOOOAU. In this simulation the perturber just 
grazes the edge of the disc (frame 1). Disc material flows 
from the primary onto the secondary via a spiral arm or 
bridge (frame 2), while a second, opposing arm is formed in 
the disc (frame 3) . As the two stars move away, the accretion 
arm between them breaks off, but the disc around the pri- 
mary star remains highly perturbed (frame 4). This m = 2 
perturbation leads to the fragmentation of the disc around 
the primary, and the formation of a binary companion to 
the primary (frame 5). This newly-formed binary system is 
surrounded by some circumbinary material in which it soon 
clears a gap (frame 6) . By the end of the encounter, the per- 
turbing star has accreted a disc of radius 475AU and mass 
O.OSMq. The primary star has a small circumstellar disc, 
of radius 60AU and mass O.IMq. Its binary companion has 
mass 0.04Mq and no circumstellar material. The two are in 
an orbit with eccentricity 0.1 and periastron 90 AU. They are 
surrounded by a circumbinary disc of mass O.ISMq and ra- 
dius 350 AU. The binary system and the perturbing star are 
in an orbit that has a periastron of 970AU and eccentricity 
0.89. 



4.2 run dsl7 : retrograde coplanar encounter 

Encounters with (j) — tt are retrograde coplanar encounters, 
in which the spins of the disc and the orbit are anti-parallel. 
The results of run dsl7 are shown in Fig. ^. This encounter 
has periastron 500AU. The perturber approaches from the 
left-hand side of the disc and orbits anti-clockwise about it. 
As the perturber passes through the disc, a trailing shock 
forms (frame 1), and the material in the shock begins to 
fragment. Because of the exchange of momentum between 
the perturber and the shocked material, the shocked ma- 
terial (including the fragments) is no longer centrifugally 
supported relative to the primary star, and it falls into the 
region around the primary star (frame 2). This causes the 
material around the primary star to spin-up and become ro- 
tationally unstable. As a consequence, it repeatedly throws 
off spiral arms, until eventually one of the arms becomes 
self-gravitating and condenses out to produce a binary com- 
panion to the primary star. By the end of the simulation, 
a large fraction of the original extended circumstellar disc 
around the primary star has been converted into a much 
more compact circumstellar disc, having mass 0.26Mq and 
radius 50AU. The primary star's binary companion has mass 
O.IOMq, and the two are on an orbit with eccentricity 0.07 
and periastron 90AU. The perturbing star has accreted a 
disc of mass O.OSMq and radius 200AU. It is bound to the bi- 
nary, with an eccentricity of 0.50 and a periastron of 320AU. 
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Figure 1. Simulation ds02 : 4> = 0, r^eri = lOOOAU. 2000 x 1500 AU region. Contour levels at [0.3, 1 (frame 2), 3, 9 (frame 4), 30, 300] 
g cm~^. The perturber approaches from the right-hand side. It a<;cretes from the disc via a spiral arm, leading to the fragmentation of 
the disc and the formation of a new protostar. The end state is a hierarchical triple system. Frame 6 shows the binary, which has cleared 
a gap in the circumbinary material. 
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Figure 2. Simulation dsl7 : 4> = w, Vperi = 500AU. 2000 X 1500 AU region. Contour levels at [0.3, 3, 30, 300] g cin"^. The perturber 
approaches from the left-hand side. As it passes through the disc, a trailing shock forms behind it. This dissipates energy, and the 
shocked material begins to fragment and falls onto the primary star. This causes the material around the star to become unstable and 
it repeatedly throws off spiral arms until it eventually fragments to form a binary companion to the primary star. The end state is a 
hierarchical triple system. ^ ^^^^ ^^^^ ^^^^^ ^^^^^ 
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Figure 3. Simulation dsl7 



: 500AU. Positions of SPH particles. The left-hand figure shows the triple system at the end 



of the simulation. The right-hand frame shows a close-up of the primary and its companion. The primary and its disc consist of 1010 
particles, whilst its companion consists of 380 particles. 



Figure g shows particle plots of the triple system at the end 
of the simulation, and of the primary star and its binary 
companion. Note that the compact discs shown here are not 
resolved. Because of gravity softening and shear viscosity, 
the particles have formed a pseudo-lattice. 



5 RESULTS 

The energy of an orbit is given by 

GMiM2(l-e) 
2^ ' ^> 

where Mi and M2 are the masses of the two bodies, and 
e and rperi the eccentricity and periastron of the orbit 
respectively. The simulations presented here have initially 
parabolic orbits, i.e. eccentricity e — 1 and energy _E = 0. A 
reduction in the eccentricity of the orbit corresponds to the 
two stars becoming bound to each other, whilst an increase 
in the eccentricity means that the two stars have become 
unbound from each other. 

Non-coplanar encounters leave the disc tilted with re- 
spect to the orbit, in agreement with previous results (Heller 
1993, Larwood et al. 1996). Penetrating encounters can twist 
the disc through up to 30°, while even for Vperi ~ "irdisc, the 
outer regions of the disc end up tilted by several degrees. 
Such a tilt can have a large effect upon the observed spec- 
trum of the disc, due to the reprocessing of light from the 
central star (Terquem & Bertout 1993, 1996). 

5.1 Disc truncation 

The prograde coplanar encounters (ds01-ds04) should be 
most damaging to the primary disc (cf. Toomre & Toomre 
1972), and one might expect them also to have the largest 
effect on the orbit. Conversely one might espect retro- 
grade coplanar encounters to have the least effect. Clarke 



& Pringle (1993, hereafter CP93) find that truncation of 
the primary disc is indeed greatest for prograde coplanar 
encounters. However, Hall, Clarke & Pringle (1996, here- 
after HCP96) find that it is retrograde coplanar encounters 
which remove most energy from the binary orbit. For pro- 
grade coplanar encounters with rperi > rdisc, a corotation 
resonance can actually lead to the transfer of energy to the 
binary orbit, causing the stars to become unbound from each 
other. 

Figure ^ shows the fraction of the disc mass that ends 
up bound to the primary star, bound to the perturber, and 
unbound, for the coplanar encounters we have simulated, 
whilst the first three frames of Fig. ^show the corresponding 
figures for the non-coplanar encounters. Frame 4 of Fig. H 
shows the size of the disc remaining around the primary 
star. Where the encounter has resulted in the formation of 
a binary surrounded by a circumbinary disc, the size of the 
circumbinary disc is used. 

Our results are, for the most part, in good general agree- 
ment with those of CP93, with the prograde coplanar en- 
counters being the most destructive to the disc, and the 
retrograde coplanar encounters the least destructive. It can 
be seen from Fig. ^ that the coplanar prograde encounters 
always lead to more mass being captured by the perturber, 
and more mass becoming unbound from the system - than 
the corresponding coplanar retrograde encounters. For en- 
counters with a prograde element and rperi > rdisc, i.e. pe- 
riastron greater than lOOOAU, there is a marked reduction 
in the disruptiveness of the encounters, so that for these 
non-penetrating encounters, those with <^ = 7r/4 actually 
truncate the disc less than the corresponding orthogonal en- 
counters. This can be attributed to the corotation resonance 
identified by HCP96. 

The two simulations that are not in agreement with 
CP93 are dsl7 and dsl8, which are retrograde coplanar en- 
counters with periastron 500 and lOOOAU respectively, i.e. 
rperi = O.brdisc and rperi = rdisc- These encounters are 
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= = 71 





0.5 1 1.5 3 0.5 1 1.5 2 



Figure 4. The fraction of the disc mass that ends up bound to the primary star, bound to the perturber and unbound, for coplanar 
prograde {<f> = 0) and coplanar retrograde (c/) = tt) encounters. 



the most destructive of all, in terms of the size of the disc 
remaining after the encounter. The reason that these en- 
counters are so destructive is that as the perturber passes 
through the disc, it generates a trailing shock. This shock 
dissipates enormous amounts of energy, and thereby triggers 
fragmentation of the disc, so that both dsl7 and dsl8 end up 
as multiple systems. The code used by CP93 did not model 
the self-gravity of the disc or the hydrodynamic forces, and 
so was unable to model the shocks. 



5.2 Energy exchange 

Figure ^ shows the change in energy of the orbit of the pri- 
mary and perturber, for both the coplanar and non-coplanar 
encounters. The energy is plotted in code units, in which the 
total binding energy of the disc is 6.6 and the binding energy 
outside of Q.^Vdisc is 0.82. We find that for non-penetrating 
{rperi > Tdisc), nou-coplanar encounters, those with (j) = tt/A 
are the most destructive, with the orthogonal and = 37r/4 
encounters having approximately the same effect. For pen- 
etrating encounters, those with a retrograde element can 
dissipate energy via a trailing shock as they pass through 
the disc, so that when rperi = Q-^rdisc, the encounter with 
= 37r/4 removes 30% more energy from the orbit than the 
corresponding tf) — n/i encounter. 

For prograde coplanar encounters, we find, in agreement 
with the results of HCP96, that encounters with periastron 
outside the disc radius lead to the transfer of energy to the 
binary orbit, and the subsequent unbinding of the system, 
whilst the retrograde encounters remove energy from the 
orbit. For the grazing encounters, with Vpsri = rdisc, the 
prograde encounter removes marginally more energy than 
the corresponding retrograde encounter. For the penetrating 
encounters, with rperi = O.Srdisc, the retrograde encounter 
is dominated by the energy loss in the trailing shock behind 
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Figure 7. Energy transfer to the binary orbit, as a function of 
inclination, for rpf,ri/'>'disc=0-5, 1.0, 1.5, 2.0 



the perturber, and so dissipates more than five times the 
energy lost in the correspo nding prograde encounter. 

McDonald & Clarke ( 1995 ) assumed in their calcula- 
tions that the energy removed from the orbit during an en- 
counter was equal to the binding energy of the disc outside 
of periastron. From the results in Fig. m, it can be seen that 
the change in orbital energy is in fact approximately twice 
the binding energy at periastron for rperi = 0.5rdisc, for all 
except the retrograde coplanar encounter, in which much 
more energy is dissipated via the trailing shock. 
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Figure 5. Frames 1-3 show the fraction of the disc mass that ends up bound to the primary star, bound to the perturber and unbound, 
for non-coplanar encounters with <p = it /A, it /2 and 37r/4 respectively. Frame four shows the size of the disc remaining around the primary 
star after the encounter for all simulations. 



5.3 Capture rates 



cr{r) = Tvr^ + 



A-KGmr 



(2) 



Figure |^ shows the same data as presented in Fig. ^, but now 
with the change in orbital energy plotted as a function of the 
inclination of the encounter, for the four different periastra 
investigated. These data can be used to calculate the cross- 
section for the formation of binary systems by capture, and 
hence to determine if the energy dissipation during disc- 
star encounters provides an effective mechanism for binary 
formation. In order to do this, we follow the method set out 



by Heller (1995) 



The cross-section for an interaction with periastron r 
for two stars each of m ass m. that have a relative velocity i; 
at large separation, is (Binney fc Tremaine 1987) 



The cross-section for capture during an encounter of incli- 
nation i and velocity v is therefore 



(3) 



where rc{v,i) is the maximum periastron for an inclination 
of i for which capture will occur at a velocity v, that is the 
largest periastron for a given i that satisfies 



(4) 



where jj. is the reduced mass of the system, and is equal to 
m/2. The value of rc for a given v and i is calculated from the 
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Figure 8. Capture rate per unit number density for disc-star 
encounters, as a function of velocity dispersion. The solid line 
shows the rate for randomly distributed inclinations, the dashes 
show the rate assuming all encounters are prograde coplanar, and 
the dots show the rate corrected for previous disruption of discs. 



data in Fig. ^ by u sing two-dimensional spline interpolation 
([Press et al. 1992[). 



If the alignments of encounters are randomly dis- 
tributed then the mean cross-section for capture is 



i) sin i di. 



(5) 



If we assume that the stars in a cluster have a Maxwellian 
velocity distribution, with a one-dimensional velocity dis- 



persion V , and that the number density of stars is n, then 
the capture rate per star is given by 



ac{v)v dv. 



(6) 



The change in Fc caused by extrapolating the data from 
the inner set of data points at rperi = 0.5rdisc to Vperi = 0, 
compared to the case when it is assumed that encounters 
inside O.bvdisc have no effect, is approximately 10%. It is 
therefore likely that the effect of these very close encounters 
is not significant, because of the low frequency with which 
they occur, and that the results obtained by extrapolating 



to 



= are not greatly in error. 



Figure |^ shows the capture rate per unit number den- 
sity, Fc, for the two cases where encounters happen at ran- 
dom inclinations, and where all of the encounters are pro- 
grade coplanar encounters. It also shows the rate Fj, as- 



suming that an encounter with 



< rat 



completely 



dissipates th e disc. This modified ra te is calculated using 
F^ = F^/Fh ( [Clarke fc Pringie 199Ta| ), where F,, is the rate 
at which encounters with rperi < raise occur, and is given 
by 



Fh 



(y{rdisc,v)v dv. 



(7) 



In fact, as can be seen in Fig. |5| encounters with r < 
Tdiac remove only about 50% of the mass from the disc. The 
corrected rate therefore gives a lower limit to the capture 
rate, as it overestimates the rate of destruction of discs by 
encounters, while the uncorrected rate takes no account of 
the destruction, and so provides an upper limit, with the 
true value lying somewhere in between. 

In order to determine whether capture via disc-star in- 
teractions is an important binary formation mechanism, we 
calculate capture rates in the same s tellar environments con- 
sidered by Clarke & Pringie (1991a). The rates are shown in 
Table ^ Even in the central O.lpc of the Trapezium cluster, 
under the most extreme conditions, capture can only give a 
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Region 


n/pc ^ 


V/km 


r^/Myr-i 


Trapezium (centre) 


10* 


1.5 


0.1 


Trapezium 


2 X 10^ 


1.5 


0.02 


Open cluster 


102 


1 


3 X 10-3 


CMC 


12 


2 


6 X 10-4 



Table 3. Capture rates for star-disc encounters in different star- 
forming environments. 
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Figure 9. Mean probability of capture during a disc-star en- 
counter for a randomly-inclined encounter, as a function of veloc- 
ity dispersion. 



binary fraction of 10% after IMyr. Within the low-density 
environment of a giant molecular cloud, capture gives a bi- 
narity of less than 0.1% after IMyr, compared to an observed 
binary frequen cy of about 80% for pre-main-sequence stars 
in such clouds (Mathieu 1994). Capture during disc-star in- 
teractions is therefore probably not an important binary for- 
mation mechanism in typical star-forming regions. 

The rates calculated assume that star formation is dis- 
tributed uniformly across those regions, and the low cap- 
ture rates are primarily due to the low number of stars that 



'^peri 1 disc 


</> = 






</' = ^ 


(j) = n 


0.5 














1 
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1.5 


1 
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Table 4. Number of companions formed by disc fragmentation 



undergo encounters in such distributed star formation. For 
instance, in the case of the GM C, only 0.5% of stars wi ll 
undergo an encounter in IMyr (Clarke & Pringle 1991a). 
If star formation is dynamically-triggered, then the stars 
formed will be much more tightly clustered than is the aver- 
age across the whole of the star-forming region. For instance, 
in simulations of star formation triggered by clump-clump 
collisions in a GMC (Chapman et al. 1992, Turner et al. 
1995, Whitworth et al. 1995), the majority of the stars that 
are formed undergo encounters. The star formation in such 
events is coeval, however, so that the initial encounters that 
occur are likely to be disc-disc rather than disc-star interac- 
tions. It is possible for disc-star encounters to occur if a star 
has been stripped of its disc by a previous encounter. 

If, in such dynamically-triggered star formation, almost 
all of the stars that are formed undergo encounters, then 
the probability of capture occuring during an encounter is 
simply pc = Fc/F/i, where Th is in this case the rate of en- 



counters with r„ 



< 2ra 



Figure |^ shows the probability 



of capture during an encounter as a function of the velocity 
dispersion. It can be seen that even if every protostar under- 
goes one disc-star encounter at speed ~ 1 km s"^, this will 
lead to a binary frequency of only a few percent. Capture 
during disc-star interactions can not therefore be an impor- 
tant binary formation mechanism, even in dense regions of 
triggered star formation. 

5.4 Fragmentation rates 

In several of the simulations, an encounter acts to trigger 
instabilities within the disc that lead to the fragmentation 
of the disc, and the formation of new companions to the 
primary star. Table ^ lists the number of new stars formed 
by disc fragmentation. The lower number of fragments for 
rperi = O.drdisc as compared to r-peri = rdisc is probably due 
to two effects, the first being that the closer encounter leaves 
the disc much more truncated, and so with less mass that 
can fragment, and the second being that the inner regions 
of the disc are not fully resolved, as discussed in Section 



2.1 



We calculate the rate of formation of new protostars, 
all of which end up as companions to the primary, using a 
method similar to that above. The data from the simulations 
are used to calculate the number of new stars formed, /(r, i) 
after an encounter at inclination i and periastron r, so that 
for instance an encounter that led to the formation of a 
binary companion to the primary star would have f — 1. 
Two-dimensional spline interpolation is used to calculate / 
for a general r, i, subject to the constraint that f{r, i) > 
for all r and i. We then average over inclinations to get the 
mean number of companions formed after an encounter at 
periastron r. 



F{r) 



f(r, i) sin i di. 



(8) 



The rate of encounters within a distance r to r + dr of the 
star is obtained by differentiating Eq. ^ to give 

/2 da(r, v) 



where 



20FV3 



dr 



v^dv, 



cy{r, v) — Tvr^ + 



4:TTGmr 



(9) 



(10) 



© 0000 RAS, MNRAS 000, 000-000 



12 H.M.J. Boffin et al. 




0,8 



0.6 - 



0.4 - 



0.2 



2.5 



0.5 



1.5 

km s"^ 



2.5 



Figure 10. Rate of formation of companion stars, per unit num- 
ber density for disc-star encounters, as a function of velocity dis- 
persion. The solid line shows the rate for randomly distributed 
inclinations, the dashes show the rate assuming all encounters 
are prograde coplanar, and the dots show the rate corrected for 
previous disruption of discs. 



so that the rate of formation of companions by fragmenta- 
tion is 



_„2/4y2 da{r,v) 
dr 



F{r)v drdv. 



(11) 



We also calculate the rate corrected for the destruction of 
discs by previous encounters, r'^. = r^/Fh, where Vh is the 



rate of encounters with rn 



The results are shown in Fig. |10[ In the same way as 
for the capture calculations, we estimate the rate of forma- 
tion of companions for different star-forming environments. 
These results are presented in Table |^. In the central regions 
of the Trapezium, fragmentation triggered by disc-star inter- 
actions can lead to a sizeable proportion of the stars ending 
up in multiple systems. If all encounters are prograde copla- 
nar, then a star will have on average 1.1 companions after 
IMyr. It is important to note, however, that in a region like 
the Trapezium, the evolution of the whole system will be 
dominated by the presence of a few massive OB stars. The 
intense radiation from these OB stars is likely to strip the 
discs from other stars in the vicinity on a very short time- 
scale. 

For lower-density regions such as a standard GMC, the 
rate of formation of binary systems by this mechanism is 
lower, due to the smaller rate of encounters. In the centre of 
an open cluster, random and prograde encounters can lead 
to 10% and 40% of stars having companions respectively af- 
ter lOMyr, which is believed to be an upper limit on the 
lifetime of protostellar discs, with lifetimes an order of mag- 
nitude smaller being more likely ( Strom, Margulis fc Strom 
1989a0 The corresponding figures for a GMC are 0.5% and 
1%. Thus in the case of an open cluster, fragmentation trig- 
gered by encounters may be an important binary forma- 



Figure 11. Mean number of companions formed in a disc-star en- 
counter, as a function of velocity dispersion. The solid line shows 
the number for randomly distributed inclinations, the dots show 
the number assuming all encounters are prograde coplanar. 



tion mechanism, but it cannot be a dominant one, while for 
GMCs it leads to a negligible binary frequency. 

Again it is necessary to note that these rates are valid 
only for distributed star formation. In dynamically-triggered 
star formation it is possible for nearly every protostar formed 
to undergo at least one encounter. The mean number of com- 
panions formed during an encounter is given by n = F//Fh 
for randomly-oriented encounters, where F^ is the rate of en- 
counters with Tperi < 'irdisc, and n — Fp/F^ if all encounters 
are coplanar prograde. These are shown in Fig. If there 
is no preferred inclination for the encounter then the mean 
number of companions after an encounter is approximately 
0.2, whilst if all protostars undergo coplanar prograde en- 
counters, the mean number of companions is 0.5. These cor- 
respond to binary fractions of 20% and 50% respectively, 
as compared to an observed pre-main-sequence binary frac - 
tion of approximately 50% (Reipurth & Zinnecker 1993). 



Thus if all protostars undergo one coplanar prograde disc- 
star encounter, the resulting binary formation by disc frag- 
mentation is consistent with the observed binary fraction. 
However, the stars formed during dynamically-triggered star 
formation will be coeval, so that both interacting stars are 
likely to possess discs. In papers II and III, we present the 
results of simulations of such disc-disc interactions. 



6 CONCLUSIONS 

We have carried out a series of simulations of encounters be- 
tween protostellar discs and passing stars. The encounters 
can lead to the removal of energy from the orbit of the en- 
counter, and can cause fragmentation of the disc either via 
spiral-arm instabilities triggered by the encounter, or, in the 
case of retrograde coplanar encounters, through the forma- 
tion of a trailing shock as the star passes through the disc. 
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Region 


n/pc ^ 


y/km 


r^/Myr-i 


Tp/Myr-l 


Trapezium (centre) 


10'' 


1.5 


0.8 


2.2 


Trapezium 


2 X 10=^ 


1.5 


0.2 


0.4 


Open cluster 


10^ 


1 


0.01 


0.04 


CMC 


12 


2 


5 X lO-"' 


1 X 10-3 



Table 5. Rate of formation of companions for star-disc encounters in different star-forming environments. is the uncorrected rate 
assuming random inclinations for encounters, Fp is the corresponding rate if all encounters are coplanar prograde. 



However, we find that the energy removed from the orbit is 
not enough to lead to the formation of binary systems by 
capture at a significant rate. In contrast, disc fragmentation 
is quite effective in creating new companions for the stars 
at the centres of discs, provided the density is high enough 
for frequent star-disc interactions, as for instance in an open 
cluster environment. 

When star formation is dynamically triggered, for in- 
stance by the collision of two clumps within a Giant Molec- 
ular Cloud, or within a shell swept up by an expanding HIT 
region or stellar wind bubble, interactions are much more 
likely to occur, and so triggered fragmentation may lead to 
a much higher binary frequency, particularly if most encoun- 
ters are prograde coplanar. Dynamically-triggered star for- 
mation is likely to be coeval, however, so that during these 
interactions both protostars will possess discs. We investi- 
gate the dynamics of such disc-disc interactions in two sub- 
sequent papers. 
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APPENDIX A: CODE TESTS 

Detailed tests of the numerical m ethod used in this paper 
were published in Watkins et al. (1996). Among these was 



the Riemann shock-tube ( god 197^ ) , which studied the abil- 
ity of the method to model bulk viscosity correctly. The re- 
sults of this test were compared to those of an identical test 
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carried out by Turner et al. ( 1995 ) using the standard arti- 
ficial viscosity of Monaghan (1989). In the implementation 
of the test (in both papers) the change in internal energy, 

was neglected. In most simulations (including those pre- 
sented in this paper), a barytropic equation of state is used, 
with the internal energy being solely a function of density, 
so that it is not necessary to calculate However, for the 
shock-tube, the equation of state is not barytropic, and ne- 
glecting this term leads to the entropy remaining constant 
across the shock, whereas in fact entropy should be gener- 
ated in the shock. Therefore we have repeated the shock- 
tube simulation including the ^ term. 

The initial conditions and method used are described 
in Watkins et al. 



( 1996 ). When using the artificial viscosity, 
the rate of change of inter nal energy, due to pressure and 
viscous forces, is given by (Monaghan 1992) 



dui 
~dt 

where 

n,, = 

and 



3 



pi 



(v. 



h-Vi 



0.01/!,2 ' 



(Al) 



(A2) 



(A3) 



where CiJ is the mean value of the sound speed at particles 
i and j, and h is the SPH smoothing length. 

For a Navier-Stokes kinematic bulk viscosity, v, the vis- 
cous contribution to the rate of change of internal energy is 
given by 



dui 



(A4) 



We find that more accurate results are obtained by using a 
mean of the contributions from particles i and ,7' in a manner 
similar to that used for the artificial viscosity (Hernquist & 
Katz 1!)89|), so that the contribution becomes 



dui 



j 



where 



(V ■ v). 



(A5) 



(A6) 



We use the /3-term from the artificial viscosity, with /3 = 1, 
and Vi = Cihi/8, so that the bulk viscosity should b e equiv- 
alen t to that of the artificial viscosity with a ~ 1 (Murray 
199e|)^he total change of internal energy in this simulation 
is therefore given by taking the sum of Eqs. Al and AS, with 
a = and (3—1. 



The results of the test are shown in Fig. 



Al 



The rar- 
efaction wave, contact discontinuity and shock can all be 
seen. The analytic solution is plotted as a dashed line. The 
introduction of the ^ term leads, to more scatter in the 

at ' 

velocity in the post-shock region, but the scatter is about 
the correct mean, and the shock and contact discontinuity 
are adequately modelled. In particular, the density changes 
are correct. 

When comparing these results with those from other 



codes, it is important to note (i) that they were obtained 
with a 3-D code (much more impressive - but essentially 
irrelevant - results obtained with 1-D codes are often pre- 
sented in this context); and (ii) that the viscosity used here 
has been fine-tuned to reproduce problems involving shear- 
flows rather than shocks. Although the results do not fit the 
analytic solution perfectly, they co mpare favourably with 
those obtained by Rasio & Shapiro (1991) (their Figure 2) 
using a comparable number of particles (~ 10''). 



APPENDIX B: HIGH RESOLUTION RUN 

We have rerun simulation ds02 with more particles to see the 
effect of increasing the resolution. The disc is now made of 
11299 particles. The results are shown in Fig. Bl. Compar- 
ing this with Fig. 1, we see that the main features in the two 
runs are similar, in particular the appearance of the spiral 
arms and the fragmentation of the disc around the primary. 
By the end of the simulation, the secondary had accreted 
0.09Mq. Because the resolution is better, however, we are 
able to resolve 2 or maybe 3 companions around the pri- 
mary. Their masses are O.117M0, 0.023Mq and O.OOTM©. 
The primary is now surrounded by a disc of only 0.0325Mq. 
Thus, due to the higher resolution, mass which was still in a 
disc around the primary in the low resolution run, has been 
accreted onto the newly formed close companion. In addi- 
tion, two other less massive companions have formed from 
the disc. From this, we conclude that our lower resolution 
runs still capture the important features of the flow, and 
that with increased resolution, the genesis of new protostel- 
lar discs would be even more efficient. 



APPENDIX C: VERTICALLY RESOLVED DISC 

As another check of the physical validity of the fragmenta- 
tion we obtain, we have rerun simulation ds02 with a ver- 
tically resolved disc. To make such a disc, we duplicated 
the evolved, single-layer disc, seven times and superposed 
the duplicates in the vertical direction so as to ensure that 
the disc was in vertical hydrostatic equilibrium. This verti- 
cally resolved disc was then evolved in isolation to remove 
any transients. The way the disc was constructed, however, 
insured that it quickly settled into a stable state. The ver- 
tically resolved disc was composed of 13713 particles. 



The results of this simulation are shown in figure CI 
Again, we see that the main features of Fig. 1 (the morphol- 
ogy of the spiral arms, and fragmentation of the disc around 
the primary) are essentially the same. After 2.58 lO"* years, 
two protostars start to form around the primary. Two other 
protostars form 1300 years later, one of which is directly 
destroyed by the primary while the other merges with one 
formed previously. 

By the end of the simulation (i.e. after 3.07 10* years), 
the secondary had accreted 0.07Mq. The two companions 
around the primary had masses of 0.033Mq and 0.016Mq. 
The primary was surrounded by a disc of 0.13Mq. These 
values are comparable with those obtained in the standard 
ds02 simulation reported in Section 4.1. The differences are 
attributable to different seed noise in the different simula- 
tions. 
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Figure Bl. Simulation ds02 at higher resolution (11301 ptcls); <^ = 0, r-peri = lOOOAU. 2000 X 1500 AU region. Contour levels at [0.3, 
1, 3, 9, 30, 300] g cm~^. The perturber approaches from the right-hand side. The evolution is clearly similar to the lower resolution one 
(Fig. 1). The end state consists, however, of 4 (maybe 5) stars instead of 3. 
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Figure CI. Simulation ds02 using a vertically resolved disc (13715 ptcls): 4> = 0, Vp^ri = lOOOAU. 2000 X 1500 AU region. Contour 
levels at [0.3, 1, 3, 9, 30, 300] g cni^"^. The pcrturber approaches from the right-hand side. The evolution is similar to the lower resolution 
one (Fig. 1). The end state consists of 4 stars. (Note that the end state of this simulation is a little earlier than for the corresponding 
two-dimensional simulations illustrated in Figures 1 and Bl.) 
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